/*
 *  main5a.cpp
 *  
 *
 *  Created by Camilo Ruiz Méndez on 29/11/08.
 *  Copyright 2008 __MyCompanyName__. All rights reserved.
 *
 */



#include <iostream>
#include <vector>
#include <string>

#include <alloc.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex>
#include <fstream.h>

#include "../include/Parameter.h"
#include "../include/carpetas.h"
#include "../include/gridobject.h"
#include "../include/laserfield.h"
#include "../include/routinesO.h"
#include "../include/constants.h"

#define complex complex<double>

using namespace std;
using std::cout;
using std::vector;
using std::string;

Parameter<string> jobname("jobname","Run1");
Parameter<double> lambda_nm("lambda_nm",32, "Wavelength in nm");
//                                          ^^ optional comment
Parameter<vector<int> > point("point", vector<int>(3,0) );

Parameter<int>    np0("np0",4000,"Particle number");
Parameter<int>    ndg0("ndg0",8192,"Spatial grid points");
Parameter<int>    ndt0("ndt0",32768,"Temporal grid points");

//Field parameters
Parameter<double> w0("w0",0.057,"Laser frequency");
Parameter<double> E0("E0",4,"Laser Amplitude");
Parameter<int>    nw("nw",20,"Number of laser cycles");

//Density
Parameter<double> density("density",15.,"Density in critical density units");

//Parameter for the average
Parameter<int> na("na",3,"Smoother");




int main()
{

  // read configuration file
  ConfigFileParser cfg( "example.cfg" );
  
  ParameterMap pMap;
  // initialize all Parameter<T> type objects from the config file
  pMap.init( cfg );
  
  
  cout<<"Wavenumber: "<< 1/( lambda_nm * 1e-9 ) <<"\n";
  
  // string/map access:
  double x; cfg.getValueInto(x, "lambda_nm"); // ConfigFileParser style
  double y = pMap["lambda_nm"];               // ParameterMap style
  string s = pMap["jobname"];                 // kind of "chameleon"
  
  try {
    cout << "Trying to pass a Parameter<double> to an int:\n";
    int i = pMap["lambda_nm"];
  } catch ( std::runtime_error e ) {
    cerr << "   ->    " << e.what() <<"\n";
  }
  
  cout<<pMap;
  
  
  
  /*************************************************/
  //Create folders and files.
  FILE *fp0; //Fichero de prueba
  FILE *fo1, *fo2; //Ficheros de salida de datos
  int dir_limit=50; //Numero de maximo de directorios
  
  //Creating a new data directory
  
  if(CrearCarpeta(dir_limit)==0)
    {	
      printf("Has superado el limite de directorios que tu pusiste, cazurro!!");
      exit(1);
    }
  
  
  fo1=fopen("output1.txt", "wt+"); //Fichero donde se almacenan los campos.
  if(fo1==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  fo2=fopen("output2.txt", "wt+"); //Fichero donde almaceno las posiciones y las velocidades.
  if(fo2==NULL)
    {
      printf("Error opening the file.\n");
      exit(2);
    }
  
  
  /*************************************************/
  // Define variables
  /*************************************************/
  
  
  
  //int iaux,npion,npim;
  int npion=0;
  
	int np=pMap["np0"];//4000; //40000;//400;//3; //4000;              ///< Particle number
	int ndg=pMap["ndg0"];//1024*8; //1024*8; //200; //1024*8;           ///< Number of points in the spatial grid
	int ndt=pMap["ndt0"];//1024*32; //1024*32; //5000;//1024*32;      ///< Number of points in the temporal grid
  
  
	//Variables for output
	int nscleft, nscright;
	//int if1,if2,if3,if4,if5;
	int nsc;
	int idj[10000];
	
	
	//Field parameters
	//double w0=0.05; //0.057;
	//int nw=20;
	//double E0=4.;
	laser laser0(pMap["w0"],pMap["nw"],pMap["E0"]);
	
	//Calculation parameters.
	double lg=(laser0.nw+1)*laser0.lambda/4; //WATCH Fix here to make deltag=c*deltat
	double tmax=(laser0.nw+1)*2*Pi/laser0.w0;
	double deltag=lg/ndg;
	double deltat=tmax/ndt;
	
	//Material parameters
	double xgleft=lg/2.-laser0.lambda/2.;
	double xgright=xgleft+laser0.lambda/1.;
	double nc=laser0.w0*laser0.w0*mass/4./Pi/q/q;  //Volumetric density
	double n0=density*nc;//15*nc;					//Volumetric density
	
	//Limits of the material sheet
    int ileft=int(xgleft/deltag);
	int iright=int(xgright/deltag);
	
	//int na=3;
	
	double qm=-q/mass;
	double ae=0.5*qm*deltat;

	grid space;
	grid time;
	space.set_grid(ndg,deltag);
	time.set_grid(ndt,deltat);
	
	
	//Arrays over the spatial grid.
	grid_object Ro(space),Rop(space),Jx(space),Jy(space),Jz(space);
	grid_object Roold(space), Ropold(space);//,p0(space);
	grid_object IntRo(space), IPhiM(space), IPhiP(space);
	grid_object IntJx(space),IntJy(space),IntJz(space);
	grid_object IAxM(space), IAxP(space),IAyM(space),IAyP(space),IAzM(space),IAzP(space);
	grid_object Ax(space),Ay(space),Az(space),phi(space), DAx(space),DAy(space),DAz(space),Axold(space),Ayold(space),Azold(space);
	grid_object Ex(space),Ey(space),Ez(space),By(space),Bz(space);
	
	
	double xg;//,deltag;//deltat,tmax,gamma;
	
	vector<double> xp(np,0.), vpx(np,0.),vpy(np,0.),vpz(np,0.);
	
	double Epx,Epy,Epz,Bpx,Bpy,Bpz;
	
	//FFT variables
	double sp1[ndt],sp2[ndt],sp3[ndt],sp4[ndt];//,raux;
	 
	//Particle parameters
	vector<double> qp(np,0.),mp(np,0.);
	
	int iflagion=0;
	int labelx=1, labely=2, labelz=3; //labels used as arguments for the function CalculaJ
	

		
    //Printing parameters in a file
    fprintf(fo1,"PARAMETROS DEL CAMPO\nE0=%e\nw0=%e\nT0=%e\nlambda=%e\nnw=%d\n\n", laser0.E0, laser0.w0, 2*Pi/laser0.w0, laser0.lambda, laser0.nw);
   // fprintf(fo1,"PARAMETROS DEL MATERIAL\nxgleft=%e\nxfright=%e\nnc=%e\nn0=%e\nileft=%d\niright=%d\n\n", xgleft, xgright, nc, n0, ileft, iright);
   // fprintf(fo1,"\nPARAMETROS DEL CALCULO\nLg=%e\ntmax=%e\nndg=%d\nndt=%d\ndeltag=%e\ndeltat=%e\n\n", lg, tmax, ndg, ndt, deltag, deltat);
   // fprintf(fo1,"PARAMETROS DE LAS PARTICULAS\nnp=%d\nna=%d\nqm=%e\n\n", np, na, qm);
   // fprintf(fo1,"n    t      Eyreflejado    Eytransm    Eyincid    Exleft    Exright    Roleft    Roright\n", np, na, qm);
    fprintf(fo2,"n(tiempo)   j(num de particula) xp     vpx            vpy         vpz\n");

	/******************************************************/
	/* Initialize variables
	/* Calculate Densities, Fields and position at the beggining
	/******************************************************/
	

	if (iflagion==0)
	{
		 npion=np;
		for(int j=0;j<np;j++)
		{
			qp[j]=-q*n0*(xgright-xgleft)/np;
			mp[j]=mass*n0*(xgright-xgleft)/np;
		}
	}
	else
	{
		npion=0;
	}
	
	//Writing parameters
	nsc=100;
	nscleft=int(xgleft-0.2*laser0.lambda)/deltag;
	nscright=int(xgright+0.2*laser0.lambda)/deltag;
	
	for(int j=0;j<nsc;j++)
	{
		idj[j]=nscleft+int((nscright-nscleft)/(nsc-1.)*(j-1.));
	}
	

	//Initialize particles at rest.	
	if(iflagion==0)
	{
		for(int i=0;i<npion;i++)
		{
			xp[i]=(xgright-xgleft)/npion*1.*(i+1)+xgleft;
            vpx[i]=0.;
			vpy[i]=0.;
			vpz[i]=0.;
			double gamma=1./sqrt(1.-(vpx[i]*vpx[i]+vpy[i]*vpy[i]+vpz[i]*vpz[i])/lightC_au/lightC_au);
			vpx[i]=vpx[i]*gamma;
			vpy[i]=vpy[i]*gamma;
			vpz[i]=vpz[i]*gamma;
		}
	}//end if
	
	//Calculate the velocities;
 	CalculaRo(xp,Ro,qp, deltag, na, fp0);
	CalculaJ(labelx,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labelx is 1 and means 'x'
	CalculaJ(labely,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labely is 2 and means 'y'
	CalculaJ(labelz,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labelz is 3 and means 'z'
	
	//Density of positive charge
	if(iflagion==0)
	{
		for(int i=0;i<ndg;i++)
			Rop.v[i]=-Ro.v[i];
	}
	else
	{
		for(int i=0;i<ndg;i++)
		{
			Rop.v[i]=0.;
		}		
	}//end if
	

	
	/******************************************************/
	/* Central part of the program
	/* Time and particle loop 
	/******************************************************/
	int pacer=0;
	for(int n=0;n<ndt;n++)
	{
		if(n%100==0) printf("n=%d\n", n);
        pacer+=1;
		double t=n*deltat;
		for(int j=0;j<npion;j++)
		{
			
            interpola(xp[j],Ex,Ey,Ez,By,Bz,Epx,Epy,Epz,Bpx,Bpy,Bpz,t,ndg,deltag,na, fp0);
			
			//Boris algorithm to solve the equation of motion
            ParticleMover(xp[j],vpx[j],vpy[j],vpz[j],Epx,Epy,Epz,Bpx,Bpy,Bpz,deltat,ae);
			
			//Printing positions and velocities
			if(j%40==0 && n%(ndt/1024)==0)
			          fprintf(fo2,"%d\t%d\t%e\t%e\t%e\t%e\n",n, j, xp[j], vpx[j], vpy[j], vpz[j]);
		
		}//End of the Boris algorithm to solve the equation of motion
		
        CalculaRo(xp,Ro,qp, deltag, na, fp0);
        CalculaJ(labelx,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labelx is 1 and means 'x'
	    CalculaJ(labely,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labely is 2 and means 'y'
        CalculaJ(labelz,xp,vpx,vpy,vpz,Jx,Jy,Jz,qp, deltat, deltag, na, fp0, fp0); //labelz is 3 and means 'z'
		
		
		
		/******************************************************/
		/* Calculate Phi, Ax, Ay and Az
		/* By solving the wave equation 
		/******************************************************/


		for(int i=1;i<ndg-1;i++)
		{
			//Time integrals of Ro and J
			
			IntRo.v[i]=IntRo.v[i]+(Ro.v[i]+Rop.v[i]+Roold.v[i]+Ropold.v[i])*0.5*deltat;  //Ros and j are defined in diff time positions
			
			IntJx.v[i]=IntJx.v[i]+Jx.v[i]*deltat;
			IntJy.v[i]=IntJy.v[i]+Jy.v[i]*deltat;
			IntJz.v[i]=IntJz.v[i]+Jz.v[i]*deltat;
			
			phi.v[i]=2.*Pi*lightC_au*(IPhiM.v[i-1]+IPhiP.v[i+1]+deltag*IntRo.v[i]);
			
			Ax.v[i]=2.*Pi*(IAxM.v[i-1]+IAxP.v[i+1]+deltag*IntJx.v[i] );
			Ay.v[i]=2.*Pi*(IAyM.v[i-1]+IAyP.v[i+1]+deltag*IntJy.v[i] );
			Az.v[i]=2.*Pi*(IAzM.v[i-1]+IAzP.v[i+1]+deltag*IntJz.v[i] );
		}
		
		IntRo.v[0]=IntRo.v[0]+(Ro.v[0]+Rop.v[0]+Roold.v[0]+Ropold.v[0])*0.5*deltat;
		
		IntJx.v[0]=IntJx.v[0]+Jx.v[0]*deltat;
		IntJy.v[0]=IntJy.v[0]+Jy.v[0]*deltat;
		IntJz.v[0]=IntJz.v[0]+Jz.v[0]*deltat;
		
		phi.v[0]=2.*Pi*lightC_au*(IPhiP.v[1]+deltag*IntRo.v[0]);
		Ax.v[0]=2.*Pi*(IAxP.v[1]+deltag*IntJx.v[0]);
		Ay.v[0]=2.*Pi*(IAyP.v[1]+deltag*IntJy.v[0]);
		Az.v[0]=2.*Pi*(IAzP.v[1]+deltag*IntJz.v[0]);
		
		IntRo.v[ndg-1]=IntRo.v[ndg-1]+(Ro.v[ndg-1]+Rop.v[ndg-1]+Roold.v[ndg-1]+Ropold.v[ndg-1])*0.5*deltat;
        IntJx.v[ndg-1]=IntJx.v[ndg-1]+Jx.v[ndg-1]*deltat;
		IntJy.v[ndg-1]=IntJy.v[ndg-1]+Jy.v[ndg-1]*deltat;
		IntJz.v[ndg-1]=IntJz.v[ndg-1]+Jz.v[ndg-1]*deltat;
		
		phi.v[ndg-1]=2.*Pi*lightC_au*(IPhiP.v[ndg-2]+deltag*IntRo.v[ndg-1]);
		Ax.v[ndg-1]=2.*Pi*(IAxP.v[ndg-2]+deltag*IntJx.v[ndg-1]);
		Ay.v[ndg-1]=2.*Pi*(IAyP.v[ndg-2]+deltag*IntJy.v[ndg-1]);
		Az.v[ndg-1]=2.*Pi*(IAzP.v[ndg-2]+deltag*IntJz.v[ndg-1]);
		
		for (int i=0;i<ndg;i++)
		{
			DAx.v[i]=(Ax.v[i]-Axold.v[i]);
			DAy.v[i]=(Ay.v[i]-Ayold.v[i]);
			DAz.v[i]=(Az.v[i]-Azold.v[i]);			
		}
		
		//Spatial integrals.
		for (int i=0;i<ndg-1;i++)
		{
			IPhiP.v[i]=IPhiP.v[i+1]+IntRo.v[i]*deltag;
			IAxP.v[i]=IAxP.v[i+1]+IntJx.v[i]*deltag;
			IAyP.v[i]=IAyP.v[i+1]+IntJy.v[i]*deltag;
			IAzP.v[i]=IAzP.v[i+1]+IntJz.v[i]*deltag;
		}
		
		IPhiP.v[ndg-1]=IntRo.v[ndg-1]*deltag;
		IAxP.v[ndg-1]=IntJx.v[ndg-1]*deltag;
		IAyP.v[ndg-1]=IntJy.v[ndg-1]*deltag;
		IAzP.v[ndg-1]=IntJz.v[ndg-1]*deltag;
		
		for(int i=ndg-1;i>=1;i--)
		{
			IPhiM.v[i]=IPhiM.v[i-1]+IntRo.v[i]*deltag;
			IAxM.v[i]=IAxM.v[i-1]+IntJx.v[i]*deltag;
			IAyM.v[i]=IAyM.v[i-1]+IntJy.v[i]*deltag;
			IAzM.v[i]=IAzM.v[i-1]+IntJz.v[i]*deltag;
		}

		IPhiM.v[0]=IntRo.v[0]*deltag;
		IAxM.v[0]=IntJx.v[0]*deltag;
		IAyM.v[0]=IntJy.v[0]*deltag;
		IAzM.v[0]=IntJz.v[0]*deltag;
		
		for(int i=0;i<ndg;i++)
		{
			Axold.v[i]=Ax.v[i];
			Ayold.v[i]=Ay.v[i];
			Azold.v[i]=Az.v[i];

			Roold.v[i]=Ro.v[i];
			Ropold.v[i]=Rop.v[i];
		}
		
		/******************************************************/
		/* End of the
		 /* Calculation Phi, Ax, Ay and Az
		 /* By solving the wave equation 
		 /******************************************************/
		

		
        CalculaEx(Ex,phi,DAx,t+deltat,deltat, fp0);
		CalculaEy(Ey,DAy,t+deltat,deltat,xg,xgleft,laser0);//wp,kp,E0,k0,w0);
		CalculaEz(Ez,DAz,t+deltat,deltat);		
		CalculaBy(By,Az,t+deltat,deltat);
		CalculaBz(Bz,Ay,t+deltat,deltat,xg,xgleft,laser0);//wp,kp,E0,k0,w0);			
	
		
		
		if(iflagion!=0)
		{
			//Ionization();
			//Calculap0();
		}
		
		
		//Calculate the transmitted and reflected field.
		sp1[n]=Ey.v[int(((xgleft-laser0.lambda)/deltag)-1)]-E0y( (xgleft-laser0.lambda)-1,t+deltat,xgleft,laser0);//wp,kp,E0,k0,w0);
		sp2[n]=Ey.v[int(((xgright+laser0.lambda)/deltag)-1)];
		
		//Printing in fo1, output1
        if(n%(ndt/1024)==0)
        {
                   fprintf(fo1, "%d\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n", n, t/2/Pi*laser0.w0, sp1[n], sp2[n],E0y(xgleft-laser0.lambda,t+deltat,xgleft,laser0), Ex.v[ileft], Ex.v[iright], Ro.v  
		[ileft], Ro.v[iright]);
    }
	}//end time loop
	  
	
	fclose(fo1);
    fclose(fo2);
	
	
}

void f1(grid a)
{
	cout << a.dx1;
}
